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^^ with other ecosystems, its species composition is resihent against small disturbances but strong perturba- 

I tions such as antibiotics can affect the consortium dramatically. Antibiotic cessation does not necessarily 

restore pre-treatment conditions and disturbed microbiota are often susceptible to pathogen invasion. 

[~~ ^ Here we propose a mathematical model to explain how antibiotic-mediated switches in the microbiota 

^^ composition can result from simple social interactions between antibiotic-tolerant and antibiotic-sensitive 

^"^ bacterial groups. We build a two-species (e.g. two functional-groups) model and identify regions of dom- 

ination by antibiotic-sensitive or antibiotic-tolerant bacteria, as well as a region of multistability where 
domination by either group is possible. Using a new framework that we derived from statistical physics, 

,.0 we calculate the duration of each microbiota composition state. This is shown to depend on the balance 

between random fluctuations in the bacterial densities and the strength of microbial interactions. The 
singular value decomposition of recent metagenomic data confirms our assumption of grouping microbes 
as antibiotic-tolerant or antibiotic-sensitive in response to a single antibiotic. Our methodology can be 
extended to multiple bacterial groups and thus it provides an ecological formalism to help interpret the 

/-Y-^ present surge in microbiome data. 

00 

04 Author Summary 

/^ Recent applications of metagenomics have lead to a flood of novel studies and a renewed interest in the role 

OsJ of the gut microbiota in human health. We can now envision a time in the near future where analysis of 

7—i microbiota composition can be used for diagnostics and the rational design of new therapeutics. However, 

^ most studies to date are exploratory and heavily data-driven, and therefore lack of mechanistic insights 

•^ on the ecology governing these complex microbial ecosystems. In this study we propose a new model 

<^\ grounded on ecological and physical principles to explain intestinal microbiota dynamics in response 

S to antibiotic treatment. Our model explains a hysteresis effect that results from the social interaction 

between two microbial groups, antibiotic-tolerant and antibiotic-sensitive bacteria, as well as the recovery 

allowed by stochastic fluctuations. We use singular value decomposition for the analysis of temporal 

metagenomic data, which supports the representation of the microbiota according to two main microbial 

groups. Our framework explains why microbiota composition can be difhcult to recover after antibiotic 

treatment, thus solving a long-standing puzzle in microbiota biology with profound implications for 

human health. It therefore forms a conceptual bridge between experiments and theoretical works towards 

a mechanistic understanding of the gut microbiota. 
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Introduction 

Recent advances in metagenomics provide an unprecedented opportunity to investigate the intestinal 
microbiota and its role in human health and disease [ll[2]. The analysis of microflora composition has 
a great potential in diagnostics [3] and may lead to the rational design of new therapeutics that restore 
healthy microbial balance in patients pjjG]. Before the clinical translation of human microbiome biol- 
ogy is possible, we must seek to thoroughly understand the ecological processes governing microbiota 
composition dynamics and function. 

The gastro-intestinal microbiota is a highly diverse bacterial community that performs an important 
digestive function and, at the same time, provides resistance against colonization by entero-pathogenic 
bacteria [7|]9] . Commensal bacteria resist pathogens thanks to resources competition l][8 , growth inhibi- 



tion due to short-chain fatty acid production 10 , killing with bacteriocins 11 12 and immune responses 



stimulation 13 14 . However, external challenges such as antibiotic therapies can harm the microbiota 



stability and make the host susceptible to pathogen colonization 15-20 . 

Despite its importance to human health, the basic ecology of the intestinal microbiota remains unclear. 
A recent large-scale cross-sectional study proposed that the intestinal microbiota variation in humans is 
stratified and fits into distinct enterotypes, which may determine how individuals respond to diet or drug 



intake 21. Although there is an ongoing debate over the existence of discrete microbiome enterotypes 22 



they could be explained by ecological theory as different states of an ecosystem 23 . Ecological theory 
can also explain how external factors, such as antibiotics, may lead to strong shifts in the microbial 
composition. A recent study that analyzed healthy adults undergoing consecutive administrations of 
the antibiotic ciprofloxacin, showed that the gut microbiota changes dramatically by losing key species 
and can take weeks to recover 1241. Longitudinal studies, such as this one, suggest that many microbial 



groups can have large and seemingly random density variations in the time-scale of weeks 25 26 . The 



observation of multiple microbial states and the high temporal variability highlight the need for ecological 



frameworks that account for basic microbial interactions, as well as random fluctuations i27-29 . 

Here we propose a possible model to study how the intestinal microbiota responds to treatment with a 
single antibiotic. Our model expands on established ecological models and uses a minimal representation 
with two microbial groups [30| representing the antibiotic-sensitive and antibiotic-tolerant bacteria in the 
enteric consortium (Fig. 1). We propose a mechanism of direct interaction between the two bacterial 
groups that explains how domination by antibiotic-tolerants can persist even after antibiotic cessation. 
We then develop a new efficient framework that deals with non-conservative multi-stable field of forces 
and describes the role played by the noise in the process of recovery. We finally support our model 
by analyzing the temporal patterns of metagenomic data from the longitudinal study of Dethlefsen and 
Relman ,24, . We show that the dynamics of microbiota can be qualitatively captured by our model and 
that the two-group representation is suitable for microbiota challenged by a single antibiotic. Our model 
can be extended to include multiple bacterial groups, which is necessary for a more general description 
of intestinal microbiota dynamics in response to multiple perturbations. 

Results 

Mathematical model 

We model the microbiota as a homogeneous system where we neglect spatial variation of antibiotic- 
sensitive (s) and antibiotic-tolerant (t) bacterial densities. Their evolution is determined by growth 
on a substrate and death due to natural mortality, antibiotic killing and social pressure. With these 
assumptions, we introduce, as a mathematical model, two coupled stochastic differential equations for 
the density of sensitives and tolerants (p^ and pt) normalized with respect to the maximum achievable 
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Figure 1. The two-group model of the intestinal microbiota with antibiotic-sensitive and 
antibiotic-tolerant bacteria. Antibiotic sensitives can inhibit the growth of tolerants and both groups 
compete for the same growth substrate. Model parameters tg and e^ represent the antibiotic sensitivity 
of sensitive and tolerant bacteria (where eg < e*), 'm,s and mt represent their affinities to substrate and 
"0 represents the inhibition of tolerants by sensitives. 



microbial density: 



dps _ Ps 
dt Ps + fpt 

dpt _ fpt 
dt Ps + fpt 



eps + Ut) ^ Fs{p) + Ut) 
ijpsPt - Pt + Ct{t) = Ft{p) + ^t{t) 



(1) 
(2) 



In the physics literature these types of equations represent stochastic motion in a non-conservative force 
field F. The first terms in F correspond to the saturation growth terms representing the indirect compe- 
tition for substrate and depend on /, which is the ratio of the maximum specific growth rates between 
the two groups. If / > 1 tolerants grow better than sensitives on the available substrate and the reverse is 
true for / < 1. They effectively describe a microbial system with a growth substrate modeled as a Monod 



kinetic 31 in the limit of quasi-steady state approximation for substrate and complete consumption from 

Both groups die with different susceptibility in response to the 
e defines the ratio of the combined effect 

While 



the microbes (see Methods for details). 

antibiotic treatment, which is assumed to be at steady-state. 

of antibiotic killing and natural mortality rates between the two groups (see Methods for details), 
the system can be studied in its full generality for different choices of e, we consider the case of e > 1 
because it represents the more relevant case where sensitives are more susceptible to die than tolerants in 
the presence of the antibiotic. A possible e{t) that mimics the antibiotic treatments is a pulse function. 
With this, we are able to reproduce realistic patterns of relative raise (fall) and fall (raise) of sensitives 
(tolerants) due to antibiotic treatment as we show in Fig. S4 in the Supplementary Information (SI) Text. 
Additionally, we introduce the social interaction term between the two groups, tpptPs, to implement com- 
petitive growth inhibition 



13 32 . In particular, we are interested in the case where the sensitives can 



inhibit the growth of the tolerants (V' > 0) , which typically occurs through bacteriocin production 33 



Finally we add a stochastic term ^ that models the effect of random fluctuations (noise) , such as random 
microbial exposure, which we assume to be additive and Gaussian. The analysis can be generalized to 
other forms of noise such as multiplicative and coloured. 



Antibiotic therapy produces multistability and hysteresis 

Wc first analyzed the model in the limit of zero noise, ^ = 0. In this case, we were interested in studying 
the steady state solutions that correspond to the fixed-points of equations (T|2 ) and are obtained imposing 



F = 0. We found three qualitatively-distinct biologically meaningful states corresponding to sensitive 
domination, tolerant domination and sensitive-tolerant coexistence (see SI Text). We evaluated the 
stability of each fixed point (see SI Text) and identified three regions within the parameter space (Fig. 2A). 
In the first region the effect of antibiotics on sensitive bacteria is very low (fe < 1) and domination by 
sensitives is the only stable state (sensitives monostability) . In the second region the effect of the antibiotic 
on sensitives is stronger than their inhibition over tolerants {fe > 1 -|- ip/e) and the only stable state is 
domination by tolerants (tolerants monostability). Finally, in the third region (1 < /e < 1 -|- V'/^) both 
sensitive and tolerant dominations are possible and stable, while the third coexistence fixed point is 
unstable (bistability) (see SI Text). This simple analysis shows that multistability can occur in a gut 
microbiota challenged by an antibiotic where one group directly inhibits the other (i.e. through the ^p 
term). Furthermore, it suggests that multistability is a general phenomenon since it requires only that 
antibiotic-sensitive and antibiotic-tolerant bacteria have similar affinities to nutrients. This is a realistic 



scenario because tolerants, such as vancomycin resistant Enterococcus 18 , are often closely related to 
other commensal but antibiotic-sensitive strains and therefore should have similar affinity to nutrients. 
Finally, the solution of equations (II]) and pi) reveals that hysteresis is present for values of fe and/^/i 
leading to multistability (Fig. 2B). Similarly to magnetic tapes, such as cassette or video tapes, which 
remain magnetized even after the external magnetic field is removed (i.e. stopping the recording), a 
transient dose of antibiotics can cause a microbiota switch that persists for long time even after antibiotic 
cessation. 

Noise alters stability points 

The previous analysis shows the existence of multistability in the absence of noise. However, the influx 
of microbes from the environment and/or the intra-population heterogeneity are expected in realistic 
scenarios and affect the bacterial density evolution in a non-deterministic fashion. This raises the question 
of how the noise alters the deterministic stable states and their stability criteria. We assume that the 
noise is a fraction of bacteria |(t) that can be added (or removed) at each time step, but on average 
has no effect since (^) = 0. This assumption is justified by the fact that a persistent net flux of non- 
culturable bacteria from the environment is unrealizable. We also assume that this random event at time 
t is not correlated to any previous time t' , which corresponds to {^k{t)£,k'{t')) = DS{t — t')6kk', where 
D characterizes the noise amplitude and S is the Dirac delta function. We calculated the stationary 
probability of the microbiota being at a given state by solving the stationary Fokker-Planck Equation 
(FPE) ^34! corresponding to the Langevin equations (1,2): 

- V • (F P,) + ^V^P, = 0. (3) 

By numerically solving equation ([3]) as described in 35 , for increasing D, we find that for small values 
of D the most probable states coincide with the deterministic stable states given by F = (Fig. 3A). 
However, by increasing D the distribution Pg spreads and the locations of the most probable states change 
and approach each other. As a consequence, the probability of an unstable coexistence, characterized by 
Ps > and pt > 0, increases thus avoiding extinction. This intuitively justifies how recovery to a sensitive- 
dominated state within a finite time after antibiotic cessation becomes possible with the addition of the 
noise. Without noise, the complete extinction of sensitive bacteria would have prevented any possible 
recolonization of the intestine. Beyond a critical noise level (Dc) bistability is entirely lost and the 
probability distribution becomes single-peaked with both bacterial groups coexisting. The microbiota 
composition at the coexistence state can be numerically determined from the solution of Ps{ps,Pt), as 
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Figure 2. Multistability and hysteresis in a simple model of the intestinal microbiota. A: phase 
diagram showing the three possible stability regions. Antibiotic-sensitive bacteria dominate when 
/e < 1 and antibiotic-tolerant bacteria dominate when fe > 1/2 -|- ^1 -I- AfTp/2 and therefore these are 
regions of monostability. There is a region of bistability between the two regions where domination by 
either sensitives or tolerants is possible. B: schematic display of the hysteresis phenomenon explaining 
cases where antibiotic treatment produces altered microbiota (i.e. tolerants domination) that persists 
long after antibiotic cessation. C-F: mean density values obtained simulating the Langevin dynamics for 
a maximum time T = 10000 after an instantaneous change of the parameter ij: (C and D) and e (E and 
F). These averages are obtained over 1000 noise realizations. C, D and E, F show the 
antibiotic-tolerants or antibiotic-sensitives densities, respectively, as a function of the social interaction 
parameter [ij:) with fe ~ 1.21 or the antibiotic killing (e) with fij] = 0.77. 



shown in Fig. 3B and Supplementary Video SI. Further investigations based on analytical expansion 
of the Langevin equations (see Methods) show that for small random fluctuations, D <C D^, the first 
noise-induced corrections to the deterministic density are linearly dependent on D with a proportionality 
coefficient determined by the nature of the interactions (insets in Fig. 3B). These linear correction terms 
can be obtained as a function of the model parameters and, after substituting a particular set of values 
in the bistable region (/ = 1.1, e ^ 1.1 and tp — 0.4), they are (Q ) = —4.3 D for sensitives and 
(Q ) = 4.4 D for tolerants. These numbers are different from those reported in the insets of Fig. 3B. 
However the discrepancy is due to the propagation of the boundary conditions when numerically solving 
the solution of the FPE using finite elements (see SI Text). 

This has important biological implications since it suggests that extinction is prevented and, more 
importantly, that a minority of environmental microbes can settle in the gut at a rate that depends on 
the strength of their social interaction with the established microbiota. 

The introduction of random perturbation affects the stability criteria of the stable states. In particular, 
we observe that the bistability region decreases when the noise amplitude D increases (Fig. 2C-F). At 
the limit, when D > Dc the bistability is entirely lost and the only stable state is the one where both 
groups coexist. This concept was previously hypothesized but not explicitly demonstrated in a model of 



microbial symbionts in corals 30 



Noise affects the recovery time 

Our model predicts that in absence of stochastic fluctuations the recovery time is larger than any obser- 
vational time-scale so that it is impossible to revert to the conditions preceding antibiotic perturbation 
(see Fig. S4 in SI Text). In reality, data show that this time can be finite and depends on the microbiota 
composition and the degree of isolation of the individuals 18 24 36 . Thus, we aim to quantitatively 
characterize how the relative contribution of social interaction and noise level affects the computation of 
the mean residence time. 

In order to determine the relative time spent in each domination state, we compute the probability of 



residence 7ri(t) in each stable state i — 1,2, ..,N using master equations 34 . This method is more efficient 
than simulating the system time evolution by direct integration of the Langevin equations because it boils 
down to solving a deterministic second-order differential equation. Furthermore, this approach scales up 
well when the number of microbial groups increases, in contrast to the numerical solution of the FPE 
which can become prohibitive when A^ > 3. In our model, the master equations for the probability ni{t) 
of residing in the tolerant i = 1 or sensitive i = 2 domination state are: 

"^ = -V2^l7r2{t) + Vi^2Mt) (4) 

at 

where Vi^j is the transition rate from state i to j, which can be obtained in terms of the sum over all 
the state space trajectories connecting i to j. 

By solving this system of equations at steady-state, we obtain the residence probabilities tti — 

(1 -|- 7-'i_>.2/'P2-j-i)~ and 7r2 = (1 -I- 'P2->-i/'Pi^2)~ ■ After computing the transition rate Vi-fj oc e n— 
as a function of the parameters, as reported in the Methods, we determine tt2, which is our theoretical 
prediction for the mean relative residence time (^2/(^1 + ^2)) spent in the tolerant domination state (see 
Fig. 5). The theoretical predictions are in good agreement with those obtained by simulating the dy- 
namics multiple times and averaging over different realizations of the noise. A first consequence from this 
analysis is that the time needed to naturally revert from the altered state depends exponentially on the 
noise amplitude (l/D). As such, we predict that for the case of an isolated system {D ^ 0) the switching 
time is exponentially larger than any other microscopic scale and the return to a previous unperturbed 
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Figure 3. Most probable microbiota states change from bistable scenario to mono-stable coexistence 
with increasing noise. A: the bacterial density joint probability distribution determined by solving the 
Fokker-Planck equation (pi) for four different values of the environmental noise. B: the bacterial 
densities at the peaks of P{ps,pt) as a function of the noise parameter D. Red symbols are data from 
the numeric solution of the Fokker-Planck equation and the black solid lines are the exponential fit. 
Parameters used: / = 1.1, e = 1.1 and ?/> — 0.4. The insets detail the linear regime. 



state is very unlikely. On the contrary, as the level of random exposure D is increased, the time to 
recover to the pre-treated configuration decreases (see Fig. S4 in SI Text). Additionally, this method 
can be considered as a way to indirectly determine the strength of the ecological interactions between 
microbes which can be achieved by measuring the amount of time that the microbial population spends 
in one of the particular microbiota states. Therefore, it can potentially be applied to validate proposed 
models of ecological interactions by comparing residence times measured experimentally with theoretical 
predictions. 
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Figure 4. Microbiota resident time in antibiotic-tolerant domination in function of the antibiotic 
action (e) and the social interaction (-0) parameters. Blue circles show the theoretical predictions 
obtained by determining the probability of the most probable path. Red circles are obtained by 
simulating the Langevin dynamics over 10000 iterations and averaged for 1000 noise realizations. 
Higher order-corrections can be included to increase the theoretical estimation accuracy. 



Analysis of metagenomic data reveals antibiotic-tolerant and antibiotic-sensitive 
bacteria 

We now focus on the dynamics of bacteria detected in the human intestine and test the suitability of our 
two-group representation by re-analyzing the time behaviour in the recently published metagenomic data 
of Dethlefsen and Relman 



24 



The data consisted of three individuals monitored over a 10 month period 
who were subjected to two courses of the antibiotic ciprofloxacin. Since the data are noisy and complex, 
and the individual subjects' responses to the antibiotic are distinct 24 , identifying a time behaviour 



by manual screening is not a trivial task. We do it by using singular value decomposition (SVD) to 
classify each subject p phylotype-by-sample data matrix X^ into its principal components. Because of 
inter-individual variability we obtain, for each subject, the right and left eigenvectors associated to each 
eigenvalue. By ranking the phylotypes based on their correlation with the first two components we recover 



characteristic temporal patterns for each volunteer 37 38 



In all three subjects, we observe that, in spite of the individualized antibiotic effect, the two dominant 
eigenvalues or principal components together capture about 70% of the variance observed in the data 
(Fig. 5A-C). Invariably, the first component shows a decrease in correspondence to antibiotic treatment 
and reflects the behaviour of antibiotic-sensitive bacteria (green line in Fig. 5D-F). Conversely, the second 
component increases with the antibiotic treatment and represents antibiotic-tolerants (red line in Fig. 
5D-F). The observation that each subject's microbiota can be decomposed into two groups of bacteria 
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Figure 5. Analysis of microbiota response to tlie antibiotic ciprofloxacin from three subjects 24 using 



singular value decomposition identifies antibiotic-sensitive and antibiotic-tolerant bacteria. A-C: 
fraction of variance explained by the five most dominant components. D-F: plot of each sample 
component 1 (green) and 2 (red) coordinates versus sample time. G-I: sorting of the phylotypes 
log2-transformed abundance matrix based on the correlation within the two principal component. 
Above (below) the green dashed lines, we display the time series of the top 20 phylotypes strongly 
correlated (anti-correlated) with component 1 and anti-correlated (correlated) with 2 and dropping 
(increasing) during treatment, which we identify as sensitves (tolerants). Subject 3 (C,F,I) displays 
absence of sensitive bacteria for a prolonged period of about 50 days after the first antibiotic treatment. 
This confirms the fact that microbiota response to antibiotic can differ from subject to subject. 
Additionally, it also supports our model prediction of remaining locked in a tolerant-dominated state 
after antibiotic treatment cessation. 



with opposite responses to antibiotics supports the validity of the two-group approach used in our model. 
Classification of each individual's phylotypes as sensitive or tolerant can be obtained by determining 
their correlation with the two principal components (see SI Text) (information in the right-eigenarrays 
matrix from SVD). Bacteria correlated with component 1 are usually highly abundant before antibiotic 
treatment and drop strongly during treatment, often below detection. Vice-versa, bacteria correlated 
with component 2 are typically in low abundance before the antibiotic and increase with antibiotic 
administration (Fig. 5G-I). Interestingly, despite significant inter-individual differences in recovery time 
(Fig. 5G-I) and individualized response of each subject, the data show that in each individual the 
majority of bacteria are antibiotic-sensitive and only a small but significant fraction are tolerant to 
ciprofloxacin (see SI Text). The recognition of these time-patterns could be considered as a possible 
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tool to indirectly determine the susceptibility of non-culturable commensal bacteria to FDA-approved 
antimicrobial compounds. However, the presence of strains in the same phylotypes that display both 
behaviors in response to the drug may constitute a significant challenge for the success of this method. 

The time evolution of the phylotypes (Fig. 5G-I) qualitatively agrees with our theoretical prediction 
that after the antibiotic administration the system moves fast, meaning in a time smaller than any other 
observable time-scale, into a new stable state with less sensitives and more tolerants. Further, the data 
also suggest that the return to sensitive domination happens after a recovery-time scale that depends on 
the microbial composition. 

Discussion 

We present a model of inter-bacterial interactions that explains the effect of antibiotics and the counter- 
intuitive observation that an antibiotic-induced shift in microbiota composition can persist even after 
antibiotic cessation. Our analysis predicts a crucial dependence of the recovery time on the level of 
noise, as suggested by experiments with mice where the recovery depends on the exposure to mice with 
untreated microbiota il8j. The simple model here introduced is inspired by classical ecological modeling 



such as competitive Lotka-Volterra models 39|40 , but relies on mechanistic rather than phenomcnological 



assumptions, such as the logistic growth. Although more sophisticated multi-species models include 



explicit spatial structure to describe microbial consortia l33,i41-43 , our model is a first attempt to 
quantitatively analyze the interplay between microbial social interactions ('0) and stochastic fluctuations 
(D > 0) in the gut microbiota. We find that these two mechanisms are the key ingredients to reproduce 
the main features of the dynamics in response to antibiotic (sudden shifts and recovery). Our model 
can be easily generalized to include spatial variability and more complicated types of noise. Therefore 
we provide a theoretical framework to quantify microbiota resilience against disturbances, which is an 
importance feature in all ecosystems |44 . By introducing a new stochastic formulation, we were able 



to characterize composition switches within the context of state transition theory 45 46 , an important 
development over similar ecological models of microbial populations [30'. We present a new method 
to calculate the rate of switching between states that identifies the most likely trajectory between two 
stable states and their relative residence time, which can be subjected to experimental validation. Finally, 
we apply SVD to previously published metagenomic data [24|, which allows us to classify the bacteria 
of each subject in two groups according to their temporal response to a single antibiotic. The SVD 
method has been used before to find patterns in temporal high-throughput data, including transcription 



microarrays j37j and metabolomics 47 . Although our approach seems to capture well the main temporal 



microbiota patterns, we should note that the use of the Euclidean distance as a metric for microbiome 



analysis presents limitations and recent studies have proposed alternative choices 48-50 . We also opt 
for an indirect gradient analysis method 51 because we are interested in emergent patterns from the 
data regardless of the measurements of the external environmental variable (i.e. presence or absence of 



the antibiotic) 50 



We propose a mechanism of interaction between two bacterial groups to explain the lack of recovery 
observed in the experiments that can be validated in the near future. Although training the model with 
the available data sets would be of great interest, this will not be useful in practice because we need more 
statistical power to be predictive. However, we anticipate that a properly validated mathematical model 
of the intestinal microbiota will be a valuable tool to assist in the rational design of antibiotic therapies. 
For example, we predict that the rate of antibiotic dosage will play a crucial role. In order to let the 
microbiota recover from antibiotic treatment, it is better to gradually decrease antibiotic dosage at the 
first sign of average microbiota composition change, which has to be larger than the threshold community 



change represented by the day-to-day variability 26 , rather than waiting for tolerant-domination and 
then stopping antibiotic treatment. 

We show here the application of our theory to a two-bacterial group scenario because we are interested 
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in the microbiota response when chahenged with a single antibiotic. However, in more reahstic conditions 
the microbiota is subjected to different types of perturbations, which may drive it towards more alternative 
stable states. Our theory of the microbial-states switches characterization can be naturally extended to 
more than two states and consists of the solution of the linear system of equations ttP = 0, where tt is 
the array of probability of residing in each stable state and P is the matrix of transition rates among the 
states. 

The ongoing efforts to characterize the microbial consortia of the human microbiome can yield tremen- 



dous benefits to human health 52-55 . Within the next few years, we are certain to witness important 
breakthroughs, including an increase in the number of microbiomes sequenced as well as in sequencing 
depth. Yet, without the proper ecological framework these complex ecosystems will remain poorly un- 
derstood. Our study shows that, as in other complex microbial ecosystems, ecological models can be 
valuable tools to interpret the dynamics in the intestinal microbiota. 

Methods 

Full model and simplification 

The model introduced in equations [T] and [2] is derived from the more detailed model described below. We 
model the bacterial competition in a well-mixed system in the presence of antibiotic treatment by means 
of the following stochastic differential equations: 

dS rUsSps mtSpt 



dt -^ ^ ' B,{S + a) Bt{S- 

dps rUsSps 



dt S + a 
dpt _ mtSpt 

dt S + a 
dA 



--fAp,-Kp,+^,{t) 
-ipPsPt - Kpt+£,t{t) 



dt 



K{Ao - A) (5) 



where we account for two bacterial groups; the intestinal resident sensitive flora ps and an antibiotic 
tolerant one pt- Additionally, we also consider the substrate 5* and the antibiotic A densities. The 
antibiotic time evolution is simply a balance between inflow and outflow (i.e. no decay due to microbial 
degradation) where K is the system's dilution rate, which sets the characteristic microscopic time-scale, 
and Aq is the constant density of the incoming antibiotic, which can be time dependent. Similarly the 
substrate concentration, 5', results from a mass balance from influx and microbial consumption. As for 
the antibiotic, 5*0 is the constant density of the incoming nutrient (i.e. the concentration of resources 
coming from the small-intestine). The second and third terms in the right-hand side of the second 
equation in (|5| describe the amount of substrate consumed by bacterial growth assuming Monod kinetics 
where rris (rrit) is the maximum growth rate for sensitives (tolerants), a is the half-saturation constant for 
growth, which parametrizes the bacterial afflnity to the nutrient, and B^ {Bt) is the yield for growth for 
sensitives (tolerants). The last two equations describe how sensitives and tolerants grow on the substrate 
available and are diluted with the factor K. We mimic the effect of the antibiotic on the sensitives 
adding a term proportional to the sensitive density where the constant of proportionality ^A is the 
antibiotic-killing rate. We also introduce a direct inhibition term ^l^Ps^ which mimics the inhibition of 
sensitive bacteria on the tolerants (social interaction). Finally the Gaussian random variables ^g, ^t are 
the additive random patterns of exposure and represent the random microbial inflows (outflows) from 
(to) the external environment. 

It is convenient to scale the variables and set the dilution rate to unity {K — 1). Therefore, all the 
rates have to be compared with respect to the system characteristic dilution rate. Introducing S = S/ Sq , 
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Ps = Ps/(Bs5o), Pt = pt/{BtSo),A = A/Ao, 7 = (^07)/^, V' = ^/{KBsSo), m, = m./iT, m* - rut/K, 
a = a/S'o, ^s = £,s/{BsSqK) and ^^ = ^t/iBtS^K) and dropping the tilde symbols, we obtain the following 
dimensionless model: 



dS 
Itt 
dps_ 

dt 
dpt 

dt 
dA 
~dt 



1-5- ^^S- -^^S 



S- 



s- 



rUsS 
^ S + a 
_ rritS 
" S + a 

= 1-A 



Ps - "fAps - Ps+£,s 
Pt - IpPsPt - Pt+it 



(6) 



If we assume that the antibiotic is a fast variable compared to the microbial densities {ps,pt) (i-e. the 
time-scale at which the antibiotic reaches stationary state is smaller than that of the bacteria), we can 
solve for ^ = and obtain A — 1. If we also assume that the incoming substrate is all consumed in 
microbial growth, therefore maintaining the population in a stationary state with respect to the available 
resources, and that, similarly to the antibiotic, the resources equilibrate much faster than the bacterial 
densities (quasi-steady state assumption, ^ = 0), we obtain that: 



S 



1 



S + a rUsPs + mtpt 



(7) 



If we now define a new parameter e — (7 + 1) describing the relative ratio of the combination of antibiotic 
killing and natural mortality (i.e. wash-out) between sensitives and tolerants, the model reduces to the 
two variables model in p reported in equations (Ip). 



Effective potential and location of long-term states 

The introduction of random noise has the important consequence of changing the composition of the 
stable states (Fig. 3A) . In order to characterize this phenomenon, we expand the solution of the Langevin 
equations (lip) around one of the stable states obtaining the following set of equations for the variable 

C^ p Pi- 



dC 
dt 



E 



df\ 



dF, 



J'^^'^^dC.dC. 



CaU + ---+ e. 



(8) 



where to simplify the notation we drop the explicit time-dependence. We can easily recognize the first 
derivative of the force on the right-hand side as the Jacobian matrix computed in one of the minima 

gl^ — J(Pi). This equation can be solved order by order by defining the expansion C = C +C + ■ • ■ 

Pi , , 
and writing the equations for each order as: 



dO 



(0) 



dt 



dO 



(1) 



dt 



E 



E 



J.AP^)Ci"^+i. 



(0) 



(9) 
(10) 



Assuming that the initial condition at time zero is Ct(0) = 0, which can always be neglected for 
long-term behaviour, the solution of equation Q is 



e\t) 



Tdi'^p^*-*') 



Ut')- 



(11) 
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This means that the average location of the minima at zer o or der is not modified by the noise since 
(C ) oc {$,) = 0. By computing the solution of the equation (10 1 we similarly find that: 



c 



m-lJXb^^-^'\ 



V..,C,'^^\t')C^^\t')dt' 



The long-time average value of the first order correction now reads: 



/rJ^^'^w) - 2 /isii E V''-'\^^^^ (e^o cTit'))dt' 



(12) 



(13) 



The time integral can be easily computed assuming that the eigenvalues of J are negative, or at least 
their real part is, as it should be for stable fixed points; therefore we obtain that: 



In^ (C«(t)) = JE- [J"']„.^-. (d°)(oo)CrM> 



(14) 



an^ 



Thus, we find that the effect of random fluctuations is to correct the value of the stable points as if 
an external field, proportional to strength of the fiuctuations, was present. This field is equal to the 
mean square displacement at large time opportunely weighted by the inverse of the curvature of the bare 



potential around the stable points, 3{p^). The correlation can be now computed using equation (11) and 
reads: 



(enoo)Cf (oo)) =^hm ^ dt' j^ dt"Y^ [c 



,j(t-t') 



=•!(<-«") 



fia' 



{Ut')ia'{t")) 



Since {£,cr{'t')^a'{t")) — DS„^'6{t' — t") the previous equation simplifies to 



{eH^)e\^))=lunD 



/*'?[■ 



j(t-t') 



=J(t-t') 



(15) 



(16) 



which results in (C ) oc D. 



Theoretical estimate of the mean residence time 

The mean residence time in each state is proportional to the residence probability ttj (t) defined in equation 
(HI). To obtain it, we need to compute the transition rate Pi^j as a function of the model parameters as: 



i^j 



tf ^ti ^p. 



Vp Pip), 



(17) 



where U and tf are the initial and final time and Vp is the functional integral over the trajectory p{t). 
Each time trajectory p{t), solution of equations ( l|2 1, has an associated weight P{p), defined as: 



p{p)^ jv^p{i)5{i-p+np))- (18) 

By discretizing the time so that t = £t with £ = 1, . . . , Af and r the microscopic time step, we obtain 



that the Langevin equations can be written using the Ito prescription 56 



as: 



F{p'-')+e 



(19) 
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where we use the short notation p{(,t) = p^ and the initial vahie is p^ = p,;. The time discretization 



allows us to interpret the functional integral in equation ( 18 ) as: 

M 



•' 1=1 
Since the noise is Gaussian and white, its distribution now reads: 

This can be justified using the property of the delta-function J 5{t ~t')dt — 1 and its discrete time version 
Tj2i=i f{T)^ij = 1 so that /(r) = e^^ follows and 5{t — t') -> Sij/r. 

Using the properties of the delta function, and integrating out all ^ s, the continuous limit expression 



of equation (21) is 

P(p(t)) = e-^ (22) 

where S{p) — \ j^' dt' \p[t') — F(p)p has an intuitive interpretation in thermodynamics and it is related 
to the entropy production rate |57|. By using stationary-phase approximation, it turns out that in the 



computation of the rate defined in (17) only one path matters, p*, which is the most probable path. 
Higher order factors are proportional to the term AT = tj — ti 45 46] , and therefore simplify with 



the denominator in equation (21). This comes from the fact that several almost optimal paths can be 
constructed starting from p* . In the optimal path, the system stays in a stable state for a very long time, 
then it rapidly switches to the other stable state where it persists until t/ . By shifting the switching time 
one obtains sub-optimal paths that, at the leading order in D, give the same contribution of the optimal 
one and their number is directly proportional to AT. This leads to 

V.^M oc e-^^ jvp exp (-^ / '^^^^'pW^I^P(O) ■ (23) 



The functional Gaussian integral can be computed [45||46] and only provides a sub-leading correction 

_5(p*) 
to the saddle-point contribution resulting in the transition rate formula Vi^j on e o ^ which is reported 

in the Results section. 

We now need to determine the optimal path and its associated action S{p*). This path is defined as 

the one where the functional derivative of S is set to zero such that the initial and final states are fixed. 

This produces a set of second-order differential equations 

which can be solved imposing the initial conditions on Pj and p{ti). 

It is easy to verify that the downhill solution is p = F and it is associated with null action. Meanwhile, 
the ascending trajectory, which is the one leading to a non-zero action and hence gives the transition 
rate value, is not given by p = — F, as it would be for conservative field of forces. This means that in 
presence of a dissipative term the reverse optimal path from the minimum to the maximum is different 
with respect to the one connecting the maximum from the minimum of the landscape. 

As the last point, we want to show that the action associated to the optimal path can be further 
simplified by noticing that 

i?=i(|pp-|F(p)p)=0. (25) 
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We can easily prove this condition by showing that the time derivative dE/dt vanishes when equation 



(24) is satisfied and remembering that the optimal path connects two stable states where F = and 
p — 0. This property allows us to rewrite the action as: 

S{p*) = £' dt' [\p*{t')\^ - p*{t') ■ F{p*{t'))] . (26) 



We solved numerically the equation (24) using a trial-and-error approach. We varied the first- 
derivative at initial time in order to arrive as close as possible to the final point within some numerical 
precision. In principle the ideal trajectory connecting two stable points should be computed in the limit 
of p{ti) — )■ but this trajectory will take infinite time. We report three examples of most probable paths 
connecting the points i to j and reverse for a chosen set of p(ti) in Fig. S6 of the SI Text . 

Singular Value Decomposition 

We first rarefy the raw phylotypes counts matrix as in f24|. We then normalize the logarithm of the 
counts according to the following procedure: 1) we add one to all the phylotypes counts to take into 
account also for the non-detected phylotypes in each sample, 2) we log-transform the data and 3) we 
normalize the resulting matrix with respect to the samples averages. In formulae, the count associated 
to phylotype i in sample j for each subject p is 

, where fij — J2i=i ^og2{Rciw'^j + 1)/^ is the average value of the counts in each sample and N is the 
total number of phylotypes. Among all possible normalization schemes, we decide to subtract the column 
averages because we aim at identifying patterns within samples based on their correlation in bacterial 
composition. Indeed, the covariance matrix of the samples is proportional to {XP)'^XP, where (X^)'^ 
is the transpose matrix. SVD on the matrix X^ is thus equivalent to the principal component analysis 
(PCA) performed on the samples covariance matrix. 
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Supplementary Information legends 

Supplementary Information Text 

The Supplementary Information (SI) Text reports additional calculations, figures and details on: 1) 
model and relative stability analysis, 2) effect of random fluctuations and noise-induced dynamics and 3) 
Singular Value Decomposition. 

Video SI 

The video shows the stationary probability distributions Ps as a function of the sensitive and tolerant 
densities for increasing noise value D, which ranges from 10^^ to 10^^. For visualization purposes, the 
noise value associated to each movie frame is displayed as an increasing bar in the top panel. 

Video S2 

The video shows the time evolution of the two principal components for the three subjects from |24 . 
Empty circles represent untreated samples, asterisks represent samples during treatement 1 and filled 
circles represent represent samples during treatement 2. 



